Impact of random obstacles on the dynamics of a dense colloidal fluid 
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Using molecular dynamics simulations we study the slow dynamics of a colloidal fluid annealed 
within a matrix of obstacles quenched from an equilibrated colloidal fluid. We choose all particles 
to be of the same size and to interact as hard spheres, thus retaining all features of the porous 
confinement while limiting the control parameters to the packing fraction of the matrix, (f>m, and 
that of the fluid, 4>i- We conduct detailed investigations on several dynamic properties, including 
the tagged-particle and collective intermediate scattering functions, the mean-squared displacement, 
and the van Hove function. We show the confining obstacles to profoundly impact the relaxation 
pattern of various quantifiers pertinent to the fluid. Varying the type of quantifler (tagged-particle 
or collective) as well as </)m and (j)i, we unveil both discontinuous and continuous arrest scenarios. 
Furthermore, we discover subdiffusive behavior and demonstrate its close connection to the matrix 
structure. Our flndings partly confirm the various predictions of a recent extension of mode-coupling 
theory to the quenched-annealed protocol. 

PACS numbers: 64.70.pv, 82.70.Dd, 61.20.Ja, 46.65.-|-g 



I. INTRODUCTION 

Fluids that are brought into contact with a disordered 
porous matrix drastically change their physical proper- 
ties. Over the past two or three decades numerous ex- 
perimental and theoretical investigations have been ded- 
icated to studying this phenomenon |l|-|6[. The remark- 
able efforts that have been undertaken to obtain a deeper 
understanding of the underlying mechanisms are not only 
motivated by academic interest. Fluids confined in dis- 
ordered porous environments play a key role in a wide 
spectrum of applied problems ranging from technological 
applications over chemical engineering to biophysical sys- 
tems. Thus, a deeper understanding of why fluids change 
their properties under such external conditions is of great 
practical relevance. 

From a theoretical viewpoint, a quantitative descrip- 
tion of such systems is a formidable challenge for several 
reasons. For one, it is difficult to reliably and realistically 
represent the matrices found in experiments or techno- 
logical applications (such as aerogels or vycor). Having 
chosen some model for the confinement, an even more 
demanding task is to formulate a theoretical framework 
that is able to appropriately describe the interplay of con- 
nectivity and confinement of the pores pertinent to the 
matrix. Even after reducing such system to the so-called 
quenched-annealed (QA) model, where for simplicity the 
matrix is assumed to be an instantaneously-frozen con- 
figuration of an equilibrated fiuid, an appropriate the- 
ory is highly intricate. Madden and Glandt [3,3, and 
later Given and Stell pj|-KLl|. derived a formalism based 
on statistical mechanics that considers the system to be 
a peculiar mixture with the matrix being one of the com- 
ponents. Within their framework, Given and Stell de- 
rived Ornstein-Zernike-type integral equations (the so- 
called replica Ornstein-Zernike equations, ROZ) that — 
in combination with a suitable closure relation — allow to 
determine the static correlators and the thermodynamic 



properties of a QA system [12| - |16j . Their formalism is 
based on a twofold thermodynamic averaging procedure: 
one over all possible fiuid configurations for a fixed matrix 
realization, the other one over all possible matrix config- 
urations. This double-averaging procedure is adapted to 
the letter in computer simulation studies of fiuids con- 
fined in disordered porous materials; in terms of com- 
puting time this is a demanding enterprise. 

Until few years ago most theoretical and simulation 
studies focused on the static properties of QA systems 
such as structure and phase behavior. In contrast, com- 
paratively little effort had been devoted to the dynamic 
properties. The main obstacles inhibiting such investi- 
gations were, on one side, the lack of a suitable theo- 
retical framework that would have allowed to reliably 
evaluate the dynamic correlation functions of the fluid 
particles, and, on the other hand, limitations in compu- 
tational power for simulations when explicitly performing 
the above-mentioned double- averaging procedure. To the 
best of our knowledge, the first simulation studies ded- 
icated to the dynamic properties of fluids in disordered 
porous confinement were performed by Gallo and Ro- 
vere |17H19| . followed shortly after by Kim [20|. In this 
context it is worth mentioning similar investigations on 
the Lorentz gas model J2l| and on binary mixtures in 
which particles are characterized by a large disparity in 
size [22J,[23 or mass [2J], i.e., systems that can be viewed 
to be closely related to QA systems. 

Over the past years, major breakthroughs have been 
achieved in describing theoretically the dynamic proper- 
ties of fluids in confined in porous media. In particular, 
use was made of mode-coupling theory (MCT) [25l - [27| 
and of the self-consistent generalized Langevin equation 
(SCGLE) theory [25, both of which rely on inferring the 
dynamics of a system from its static structure. MCT 
has played over the last decades a central role in de- 
scribing the dynamic properties of fluids [29| , including 
in particular the glass transition. Krakoviack succeeded 



in deriving MCT-typc equations for the dynamic corre- 
lation functions, which require as an input the static cor- 
relation functions of the annealed fluid, which in turn 
are obtained from the ROZ equations. Subsequently, 
Krakoviack solved these equations for the specific case 
of a hard-sphere fluid confined in a disordered matrix of 
hard-sphere particles of equal size and evaluated single- 
particle and collective correlation functions of the sys- 
tem. Based on this knowledge, a kinetic diagram was 
traced out in the parameter space spanned by the pack- 
ing fraction of the fluid and that of the matrix, indicating 
the regions in which the system is in an arrested or in a 
nonarrested state. The resulting kinetic diagram is very 
rich and quite surprising: the regions of arrested and 
nonarrested states are separated by a line along which at 
small matrix packing fractions typc-B transitions and at 
intermediate matrix packing fractions type-A transitions 
occur. As a peculiar feature of this line, a re-entrant 
glass transition is predicted for intermediate matrix and 
small fluid packing fractions. Additionally, in the region 
of nonarrested states a continuous diffusion-localization 
transition is observed. However, we emphasize that the 
combined ROZ-I-MCT framework predicts ideal transi- 
tions, whereas in experiment or simulation transitions 
will always be smeared out. 

Throughout this contribution we employ the type- 
A/type-B terminology that is conventionally used to 
characterize the behavior of a dynamic correlator as an 
ideal transition occurs. This terminology refers to the 
decay pattern of the correlator as well as to the behavior 
of its long-time limit upon crossing a dynamic transition. 
In a "type-A" transition a correlator relaxes in a single 
step, and its long-time limit may assume arbitrarily small 
nonzero values. In a transition of "type-B," on the other 
hand, a correlator decays in a distinct two-step pattern; 
upon crossing the transition the second decay step is de- 
layed to infinity and the long-time limit of the correlator 
jumps from zero to a nonzero value. The type-B behavior 
is well-known from bulk glass formers [SOl and is usually 
attributed to a "cage effect" imposed upon fluid particles 
by their neighbors. 

This remarkably rich phase behavior has motivated 
simulation studies which — as a consequence of the con- 
siderable increase in available computation power — have 
meanwhile come within reach. Short accounts have been 
and 
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32|. In this contribution we elabo- 



rate on our investigations on the dynamic behavior of a 
fluid confined in a disordered porous matrix within the 
QA picture. For sir nplicity , and for congruence with the 
theoretical studies |25l - [27| . we consider a fluid of hard 
spheres confined in a disordered matrix formed by hard 
spheres of the same size. Investigations are based on 
event-driven molecular dynamics simulations. Within 
the two-dimensional parameter space we focus mainly on 
three representative pathways to put the theoretical pre- 
dictions to a thorough test: (i) a path at constant low 
matrix packing fraction, probing the existence of the con- 
tinuous type-B transition, (ii) a path at constant interme- 



diate matrix packing fraction, probing the discontinuous 
type-A transition and the re-entrant scenario, and (iii) 
a path at intermediate, constant fluid packing fraction, 
probing the diffusion-localization transition. Results are 
summarized in kinetic diagrams, where we consider — 
complementing our previous communication — different 
criteria to define dynamical arrest. We compare our sim- 
ulation results with the predictions of MCT, which, how- 
ever, is not straightforward due to the ideal nature of the 
latter. 

The paper is organized as follows: in Sec.|ll]we present 
the model and our simulation method, transferring the 
rather lengthy description of the system setup to Ap- 
pendix [X] Results are compiled in Sec. lIIII we start with 
the static structure of the system and then discuss various 
methods of how to define slowing down and the ensuing 
kinetic diagrams. In the subsequent three subsections 
we discuss the intermediate scattering functions and the 
mean-squared displacement along the three aforemen- 
tioned paths through parameter space. The final subsec- 
tion of Sec. mil shall deal with the van Hove correlation 
function in order to discuss trapping. In Sec. II VI we ex- 
tend the analysis of the intermediate scattering functions 
to discuss the arrest scenarios in more detail. The paper 
closes with concluding remarks alongside with an outlook 
on future work required to settle questions yet open. The 
aforementioned Appendix El is succeeded by Appendix [B] 
in which we elaborate on the issue of the effects of system 
size on the dynamical correlation functions. 



II. MODEL AND METHODS 

For the purpose of this work we opted for the 
"quenched-annealed" (QA) protocol since it is well- 
defined, well-tractable by theory, and involves only a 
small number of parameters. In this protocol both the 
porous confinement and the confined fluid are repre- 
sented by a set of particles, conventionally called "ma- 
trix" and "fluid" particles or — borrowing from the typol- 
ogy of metallurgy — as "quenched" and "annealed" com- 
ponent. In this, QA systems bear resemblance to mix- 
tures; however, the matrix particles are pinned at partic- 
ular positions and act as obstacles for the fluid particles. 
In this work, we exclusively consider systems containing 
one matrix and one fluid species; quantities pertaining to 
those species will be designated with the subscripts "m" 
and "f," respectively. 

Adopting the concept of immobile obstacle parti- 
cles, there are still numerous schemes for their arrange- 
ment [23, l32l - [35| . The QA protocol specifies the posi- 
tions of the matrix particles to be frozen from an equi- 
librium one-component fluid; the fluid particles merely 
move within this matrix and do not exert influence upon 
it [431 . Notably, the fluid particles may also dwell in 
"traps" — finite spatial domains entirely bounded by infi- 
nite potential walls formed by the confining matrix. 

Since for this work we are interested in a minimal ap- 




Figure 1: (Color online) Snapshot of a quenched-annealed 
mixture of hard spheres at 0m = 0.25 and (pf = 0.10. Dark 
(blue): matrix particles, light (gray): fluid particles. 



proach to the dynamics of a fluid moving in disordered 
confinement, we decided to exclusively employ hard- 
sphere (HS) interactions among the particles. This de- 
prives the system of any inherent energy scale, rendering 
the only relevant system attributes to be of geometrical 
nature. We further reduced the parameter space by re- 
quiring all particles involved to be of the same (monodis- 
perse) diameter a. This introduces a as the unit of length 
used throughout this study, and leaves the system to be 
governed by mere packing fractions, namely the matrix 
packing fraction 0m and the fluid packing fraction (f>f . As 
mentioned before, this particular choice of system was ex- 
amined in |25l - [27| . A snapshot of a typical HS-QA system 
at large matrix density is shown in Fig. [Ij 

In this work we study QA systems utilizing molecular 
simulations. Since we are primarily interested in dynam- 
ical features, we applied molecular dynamics (MD) sim- 
ulations to track the physical time evolution of a system. 
All MD computations were performed using the event- 
driven algorithm described in j36l |: modifications to that 
algorithm were essentially limited to fixing particles in 
space in order to adapt to the QA protocol. (Notably, 
this adjustment invalidates the conservation of momen- 
tum and angular momentum.) As is conventional, we em- 
ployed periodic (toroidal) boundary conditions and the 
minimum-image convention in order to mimic an infinite 
system. 

In addition to the usual thermodynamic averaging, 
(•••), the QA protocol imposes another averaging pro- 
cedure over disorder, • • • , so as to restore homogeneity 
and isotropy. Depending on the state point ((/)m, 0f ) and 
the number of particles in the system, we averaged our 



data over at least ten matrix configurations; for selected 
state points, we increased that number to as much as 40. 

A nontrivial task is to find an initial system configu- 
ration given some state point. In this work, we arc in- 
terested in systems in which the A'f fluid particles com- 
prise most of the space that is left accessible to them 
by the N,^ immobile matrix particles — i.e., high-density 
systems. Generating instances of the porous confinement 
is effortless, since for all state points of interest, 0m is 
situated well within the fluid regime when used for a 
one-component system of mobile particles. The problem 
rests in determining permissible initial positions for the 
A^f fluid particles. Creating an amorphous high-density 
configuration of bulk monodisperse hard-sphere particles 
already represents an interesting task, with various elab- 
orate algorithms having been developed to accomplish 
this task [33, l38[. However, none of these algorithms are 
suited for adaption to the setup of HS-QA systems. This 
unfortunate finding prompted us to devise a custom rou- 
tine which is described in detail in Appendix [K\ 

After successful setup of a system instance, the simula- 
tion time was set to t = and an attempt to equilibrate 
the system using the MD routine was made. The deci- 
sion whether or not the particular system realization was 
considered equilibrated was based on the mean-squared 
displacement 



fr2(t) = ^|rt(t)-rt(0)|') 



(1) 



where r* (t) is the location of a tagged (superscript "t" ) 
fiuid particle at time t. Note that the above definition 
is the one used in Sec. lIIII and hence involves an average 
over matrix configurations; of course this average has to 
be omitted when using Sr'^it) to characterize a single sys- 
tem realization. Also, note that the notion of a "tagged 
particle" does not imply the presence of an additional 
tracer particle. 

We defined a system realization to be in equilibrium 
if 6r^{t^) > Sri ^ 100 was fulfilled. We allowed at 
least t* ~ 30 000 time units for the system to equi- 
librate; for certain state points we extended i* to as 
much as the tenfold value. The conventional unit of 
time T = y/ma^JkBT used throughout this work is based 
upon the mass m of the fluid particles, their diameter a 
(sec above), and their temperature T [39| . 

Data were only compiled after equilibration was at- 
tempted, using the final configuration of the equilibration 
run as initial configuration of the data run. Generally, 
the final time for data collection, t^, was chosen to be 
equal to the final time of the equilibration attempt, to; in 
particular in some nonequilibrated systems we also used 
td > te to investigate particular features of the system. 
For most state points, we chose N = A^m+A^f — 1000 
to be the total particle number [431 ■ In order to retain 
reasonable statistics even close to 0f = 0, we increased A^ 
to as much as 12 000 for elevated 0m, such that A^f > 50 
for all state points. 
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Figure 2: (Color online) Fluid-fluid static structure factor 
S{k) for a series of ^f values at constant ^m. (a) (pm = 0.05, 
(b) 0m ~ 0.10. Error bars represent one standard deviation 
of the mean for different system realizations. 



Figure 3: (Color online) Fluid-fluid radial distribution func- 
tion g{r) for a series of 0f values at constant (/)m. (a) 
(pm = 0.05, (b) (jjm = 0.10. Error bars: see Fig. [2] 



III. RESULTS 



A. Static structure 



all state points indicates that no simple crystalline long- 
range order is present. 

Related arguments hold for the radial distribution 
function 



The most obvious properties to inspect in our systems 
are simple static structural properties; this investigation 
was already carried out in great detail in [401. In. the 
present context static properties are of use for a different 
reason: Since the particles were chosen to be of monodis- 
perse diameter a, crystallization might take place in the 
case of dilute matrices. Inspection of the radial distri- 
bution function g{r) and the static structure factor S{k) 
provides a means to determine whether or not this is the 
case. In fact, it is sufficient to investigate fluid-fluid corre- 
lations, gs{r) and Sff{k), since only for (/){ ^ ^m systems 
are prone to crystallization. Therefore, the index "ff" 
will henceforth be dropped. 

We first consider the static structure factor 
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(2) 



where r^j = |r* — r|, and r' denotes the location of 
fiuid particle i. Fig. [2] displays this function for mod- 
erate (0ni = 0.1) and for low {(pm = 0.05) matrix packing 
fractions, i.e., for systems particularly prone to crystal- 
lization. Clearly, all features of S{k) evolve smoothly 
upon increasing 0f, suggesting that no phase transition 
takes place. Moreover, the absence of sharp peaks for 



9(r) 




(3) 



where V is the system volume. Fig. [3] shows this func- 
tion for the same matrix packing fractions as above 
(0m € {0.1,0.05}). As for S{k), the features of ^(r) de- 
velop continuously, which again supports the notion that 
a phase transition is absent. 

Additional peaks in g{r) at specific positions would 
be a hallmark of crystallization. For fee- or bcc-like 
short-range order, the most pronounced additional peaks 
would be located at rfcc = <i>V2 or rbcc = ^VS, where 

$ = {(0m + 0f)/0max}^/^ and 0max = TT / (SV^) IS the VOl- 

ume fraction of close-packed monodisperse hard spheres. 
For instance, for (/),n = 0.05 and (j>f = 0.50, crystallites 
would cause a peak at rfcc — 1.56 or rbcc — 1.91. Clearly, 
peaks at these positions are absent in Fig.|3l On the other 
hand, in both panel (a) and (b) an additional peak at 
r ~ 1.8 emerges upon increasing 0f. However, this peak 
is a well-known feature of glass-forming systems [41| and 
thus provides further evidence for the absence of crys- 
tallization. We conclude that crystallites, albeit possibly 
temporarily present, do not play a major role in the sys- 
tems under investigation. 
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Figure 4: (Color online) Kinetic diagram based on single- 
particle properties (see text). Symbols: time ta needed for 
Fs(fc=7.0,i) to decay below F* = 0.1. Thick solid (blue) 
line: interpolation through points for which ta — 10^. Thick 
dotted (green) line: arrest line based on 5r^{t) from [3l| (ex- 
tended to low (f)f values). Thin dashed and dotted (black) 
lines: MCT transition lines pertaining to single-particle prop- 
erties from [2711. 
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system, we chose a mean-squared displacement 5r1 = 100 
and a simulation time t^, = 30 000; if the system obeyed 
(5r^(i,) > 5r1 then it was considered "nonarrested," and 
otherwise "arrested." The so-constructed kinetic dia- 
gram qualitatively confirmed the behavior of the single- 
particle dynamics of the system at hand as predicted by 
MCT [23, |31| . The arrest line determined from this cri- 
terion is indicated in Fig. 21 

In Fig. |4] we present a kinetic diagram of the system 
under investigation based on the self-intermediate scat- 
tering function 



F,{k,t) 



(pLft)p^-ic(o)) 
(pL(o)p-k(o)) 



(4) 



where p\it) = exp{ik • r*(t)} is the density of a tagged 
(superscript "t") fluid particle in k space after some 
time t has passed [50[ . The symbols in the figure indicate 
the times ts needed for Fs{k=7 .Q^t) to decay below the 
value F* = 0.1. The value k = 7.0 was chosen to be close 
to the main peak in S{k) in a typical high-density QA 
system — for instance, for 0,„ = 0.05 and (pi = 0.50 (cf. 
Fig. [2]) the first peak in S{k) is located at k ~ 7.2; for 
(/>m ~ 0.25 and 0f = 0.10 it is found at fc ~ 6.6. 

If we let i* ~ 10'^, discriminating whether for a state 
point ts > t* ("arrested") or tg < t* ("nonarrested") 
yields the thick solid (blue) line in Fig. SI Obviously, 
the latter is only slightly different from the thick dotted 
(green) line that is obtained from the 5r'^ (t) criterion de- 
scribed above. This is not unexpected since both Fs{k, t) 
and (5r^(i) are single-particle properties, and in the Gaus- 
sian limit it is even Fs{k,t) = exp|fc^Jr^(^)/6| j43 |. 
Thus, the validity of the approach chosen in [3l| is con- 
firmed. 

In order to complement the kinetic diagrams based 
upon single-particle properties, in Fig. [5] we present for 
the first time a kinetic diagram for a collective property 
of the system. For comparison with the theory developed 
in [25l - l27| we chose to operate with the connected part 
of the intermediate scattering function 



Figure 5: (Color online) Kinetic diagram based on collec- 
tive properties (see text). Symbols: time fc needed for 
Fc{k^7.0,t) to decay below F* = 0.1. Thick dash-dotted 
(red) line: interpolation through points for which fc — 10°. 
Thin dashed and solid (black) lines: MCT transition lines 
pertaining to collective properties from [271 ]. 



B. Kinetic diagrams 

Many features of the system under investigation can 
be understood by means of kinetic diagrams in which 
the "state" of a dynamic property is indicated in the 
plane spanned by 0m and (j)f. In an earlier work [3l[ we 
presented a kinetic diagram based on the mean-squared 
displacement Sr^ (t) . In order to classify the state of the 



Fe(fc,f) 



(^Pk(t)d>_k(0)) 
(5pk(0)<5p_k(0)) 



(5) 



where (5/0k(i) — Pk(i) ~ (pk(0)) is the fiuctuating part of 
Pk{t) = E, expjik • rf (i)} (cf. Sec. HV]). Very much like 
for the self-intermediate scattering function Fs{k,t), in 
this figure we display the times fc needed for Fc{k~7.0, t) 
to decay below the value F* ~ 0.1. The thick dash- 
dotted (red) line in Fig. [S] interpolates through points 
for which fc ~ 10". The shape of this iso-fc line clearly 
contradicts the MCT scenario [25|, ^M,^ which predicts 
a re-entrance regime (a "bending-back" of the collective 
glass transition line) for large values of (p^. 

Since the arrest line obtained from fg is situated at con- 
siderably larger (p-^-i than the diffusion-localization line 
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Figure 6: (Color online) Intermediate scattering functions for 
a series of <j>f values at fixed ^m ~ 0.05 and k = 7.0. (a) 
Connected correlator Fc{k,t), (b) single-particle correlator 
Fs{k,t). Error bars: see Fig. [2] 
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Figure 7: (Color online) Mean-squared displacement (a) and 
its logarithmic derivative (b) for a series of (j)i values at fixed 
4>in ~ 0.05. Error bars [only panel (a)]: see Fig. (2] 



predicted by MCT, it is reasonable to expect a similar 
shift for the glass transition line. Unfortunately, an esti- 
mate for such shift suggests that the re-entrant scenario 
may take place at combinations of </>„! and (/)f at which 
systems cannot be set up due to geometric constraints 
(cf. Appendix|X|. However, one should also bear in mind 
that the available MCT calculations are based on struc- 
ture factors obtained from integral equations |25l - [27| as 
it is still difficult to extract reliable blocked correlations 
from simulations [43, |4J1- Thus, it remains hard to dis- 
entangle the inherent deficiencies of MCT from those re- 
lated to the input structural data, especially at high ma- 
trix densities. 

Based on simple arguments, in [231 it is concluded that 
the glass transition line has to coincide with the diffusion- 
localization line in the limit of vanishing fluid density. 
The data presented in Fig. U] and [S] are ostensibly incon- 
sistent with this expectation. We remark, however, that 
our simulations do not provide for this special case since 
we decided to retain a nonvanishing number of fluid parti- 
cles in every system realization so as to always have finite 
fluid-fluid correlations. More importantly, in Sec. II VI we 
will provide evidence that using the present definition of 
the decay times tc and tg, the two quantities are not quite 
comparable to one another. 



C. Path I: constant low matrix density 



derivative 



In this subsection, as well as in Sec. lIIIDl and lllTEl we 
refine our findings from |3l| concerning Fs{k,t), Fc{k,t), 
and Sr'^{t). Additionally, we computed the logarithmic 



z{t) = 



d [log Sr^{t)] 
d[logt] 



(6) 



which facilitates the discussion of 5r^ (i) since in an as- 
sumed subdiffusive law of the form Sr^(t) ex <^ it rep- 
resents the momentary value of the power z. The first 
set of data we present for state points at various (pt at a 
common (f)^ = 0.05, which corresponds to a vertical line 
in the kinetic diagrams Fig. |4] and [5] that shall henceforth 
be referred to as "path I." 

In Fig. ini the time dependence of Fc{k, t) and Fs{k, t) 
is depicted. For both observables, the wave vector was 
chosen to be fc = 7.0, which is close to the first peak in 
the static structure factor (cf. Fig. ^. For elevated 0f 
an intermediate-time plateau can be identified in both 
Fs{k,t) and Fc{k,t); as (j)f is increased, this plateau ex- 
tends to longer and longer times but does not change in 
height. Moreover, Fs{k, t) and Fc{k, t) decay on the same 
time scale as they approach this type-B transition. This 
is a well-known phenomenon in bulk glass formers J30l |. 
where the collective glass transition drives the arrest of 
the individual particles. 

The bulklike behavior is also prevalent in Sr'^{t), as is 
evident from Fig. [71 For the lowest value of (j>f depicted, 
5r'^{t) crosses over directly from ballistic [z = 2) to dif- 
fusive {z = 1) behavior. Upon increasing (f>{ the ballistic 
range is followed by a distinct regime in which dr'^{t) 
grows only very slowly, which is reflected by a drop of 
z{t) below unity and as low as 0.15 for the highest 0f 
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Figure 8: (Color online) Intermediate scattering functions for 
a series of (j>i values at fixed ^m = 0.20 and k = 7.0. (a) 
Connected correlator Fc{k,t), (b) single-particle correlator 
Fs{k,t). Error bars: see Fig. [2] 




Figure 9: (Color online) Mean-squared displacement (a) and 
its logarithmic derivative (b) for a series of (j)i values at fixed 
4>in ~ 0.20. Error bars [only panel (a)]: see Fig. [2] 



considered. This is another manifestation of the cage 
effect (cf. Sec.|l|. For all such cpf, ensuingly diffusive be- 
havior is recovered in the long-time limit, vifith the time 
to approach z = 1 increasing enormously with (j)[. 



D. Path II: constant intermediate matrix density 

In this subsection we present data for state points 
at various ^f at a common 4>^ = 0.20, i.e., along an- 
other vertical line in the kinetic diagrams presented in 
Sec. IIIIBI In the following we will refer to this line as 
"path II." 

In Fig. [SI the time dependence of Fc{k=7.0,t) and 
Fs(fc=7.0,i) is depicted. As can easily be seen, the de- 
cay pattern of both correlators is quite unlike that in 
Sec. nil CI and Fig. [5] therein. Moreover, Fc{k,t) and 
Fs{k,t) deviate substantially from one another. Most 
strikingly, the time scales on which the correlators decay 
to their long-time value are vastly disparate, and increas- 
ingly so for larger values of 0f . At 0f — 0.22, for instance, 
the decay times read ic = 7.9x10° and tg ~ 5.4x10'^. 
This represents a difference of almost three decades in 
time, with the collective correlator decaying faster. 

Fc{k— 7.0, t) seems to approach a transition since the 
time it needs to decay to its long-time value increases 
substantially as (f>{ is increased. However, it is difficult 
to pinpoint for sure whether Fc{k,t) exhibits a transi- 
tion of type-A (single-step decay) or type-B (two-step 
decay, cf. Sec. IIIIC|) . For one, in simulations a type-A 
transition may show features that differ from the theo- 
retical predictions [2g|. Further, it is possible that an 



intermediate-time plateau would emerge upon increasing 
(j)f above 0.27 (the largest value of (j){ depicted in Fig.[8|). 
Note, however, that beyond (j)f = 0.27 geometry strongly 
inhibits to realize system instances (see Appendix [X| . 
We conclude that using the current scheme of simulation 
and analysis, no definite statement can be made about 
the type of transition in Fc {k, t) in this part of the kinetic 
diagram. 

Fs{k—7.0,t) shows an intermediate-time plateau for 
large values of (/)f , just as in Fig. [51 The times at which 
such plateau may be identified range from 10" < t < 10^ 
for (j)f = 0.15 up to 10° < t < 10-* for (f>{ = 0.27. In 
contrast to Fig.[SJ Fs{k,t) attains an additional nonzero 
long-time value that increases with (j)f. A tail can be 
identified even at the lowest fluid packing fraction con- 
sidered (0f = 0.05); whereas in this case no intermediate- 
time plateau is present. Therefore, for the sequence of 
state points at hand, two superposing decay behaviors 
exist in Fs{k,t): a type-A transition responsible for the 
long-time value, and another relaxation mechanism at 
higher 0f that leads to the intermediate-time plateau 
and is likely to be caused by a weak collective cage ef- 
fect. The type-A transition is probably connected to the 
continuous diffusion-localization transition predicted by 
MCT [131 and is discussed in more detail in Sec. IIII El 
andHnB 

Although not as obvious. Fig. [3] evidences that the 
mean-squared displacement 5r'^{t) differs from bulklike 
behavior as well. As is the case for (^m = 0.05, for 
0m = 0.20 ballistic behavior is followed by a regime for 
which z drops well below unity. Contrary to the bulklike 
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Figure 10: (Color online) Intermediate scattering functions 
for a series of (pm values at fixed (fi{ = 0.10 and k = 7.0. 
(a) Connected correlator Fc{k, t), (b) single-particle correlator 
Fs{k,t). Error bars: see Fig. (2] 
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Figure 11: (Color online) Mean-squared displacement (a) and 
its logarithmic derivative (b) for a series of (^m values at fixed 
(j>i = 0.10. Error bars [only panel (a)]: see Fig. (2] 



case, this observation holds also for very low (/)f , suggest- 
ing that caging by fluid particles is not the only mecha- 
nism effective in this regime. Most likely, at lovif 0f the 
drop in z is entirely due to the porous matrix, while at 
larger t^f both the pores and caging play a role. 

For low (pi , the range of decreasing z is followed by an 
almost straight increase to z = 1, i.e., a direct approach 
of the diffusive regime. Upon increasing 0f, a distinct 
intermediate-time plateau with z < 1 emerges in z(i), 
which corresponds to a subdiffusive regime in (5r^(t). The 
estimated value z ~ 0.78 in the subdiffusive regime is 
only weakly dependent on 0f; merely for (/)f > 0.26 the 
value of z seems to systematically decrease. This decline, 
however, may be due to insufficient equilibration of the 
respective systems. 



E. Path III: constant fluid density 

Lastly we turn to a selection of state points at vary- 
ing (/)„! for a relatively low fixed 0f = 0.10, i.e., along a 
horizontal straight line in the kinetic diagrams, perpen- 
dicular to paths I and II. This line will henceforth be 
called "path III." 

In Fig.[10]we present Fc{k^7.Q,t) and Fs{k^7.0,t) for 
these state points. The correlators differ substantially 
from each other and from their bulklike counterparts — 
even more than was the case in Sec. IIIIDI Most no- 
tably, Fs{k,t) attains a sizeable nonzero value as time 
approaches infinity, and does so more pronouncedly as (j)in 
increases. By contrast, for all 0i„ represented in Fig. [TU] 
Fc{k,t) decays strictly to zero. On the other hand the 



actual pattern of the decay to the long-time value is quite 
similar in Fc{k, t) and Fs{k, t): unlike in Fig. [51 both cor- 
relators decay in a single step, suggesting only type-A 
transitions to be involved. 

Comparison of Fig. [TT] with Fig. |9] and [7] reveals the 
functional form of Sr^ (i) to be more sensitive to (j)^ than 
to (pf. As observed in Sec. IIII Cl at low matrix densities 
the initial ballistic behavior is directly followed by diffu- 
sive behavior. For intermediate (pm < 0.25 (cf. Sec. IIII Dj) 
a region of very slowly increasing particle displacement 
{z < 1) succeeds the ballistic range. The value of the 
exponent z in this region decreases as (/),„ increases. Sub- 
sequently, diffusive behavior is recovered as z raises in an 
almost linear fashion. 

Upon further increasing ^m, at some (j}"^ a geometric 
percolation transition takes place in the space accessi- 
ble to the fluid particles ( "voids" ) , with a void stretching 
through the whole system at (j)^ < (pm- This transition is 
intimately connected with the diffusion-localization tran- 
sition predicted by MCT [23 and observed in our simu- 
lations. As a hallmark of this transition, 5r'^{t) does not 
approach diffusive behavior for (j),^ ~ 0.25 but instead re- 
mains at an approximately constant z ~ 0.5 for a window 
covering about three decades in time. Since for that ma- 
trix density z{t) ultimately increases beyond 0.5, one can 
expect the diffusion-localization transition to take place 
at slightly higher 0„i. An upper bound to that value is 
0ni = 0.2625, the smallest value presented that is larger 
than 0.25. At this packing fraction of obstacles the win- 
dow of constant z extends over roughly two decades in 
time, with z ultimately decreasing below that plateau. 
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Figure 12: (Color online) self-part of the direction-averaged 
van Hove function Ga{r,t) for (a) t/im ~ 0.05, t/if = 0.50 (b) 
0m = 0.20, (fi{ = 0.22. Curves are approximately logarithmi- 
cally spaced in time; some are highlighted for visual guidance. 
Error bars: see Fig. [2] 



Figure 13: (Color online) self-part of the direction-averaged 
van Hove function Gs{r,t) for (c) (pm — 0.25, <j>f — 0.10 (d) 
0m = 0.30, 0f = 0.05. Curves are approximately logarithmi- 
cally spaced in time; some are highlighted for visual guidance. 
Error bars: see Fig. (2] 



Judging from these observations, the subdiffusion ex- 
ponent z at the difFusion-locahzation transition might 
be sHghtly lower than 0.5. Nevertheless, the observed 
value of z is in striking agreement with MCT, which 
predicts that z = 0.5 [23|. An additional agreement of 
numbers calls for investigation: It has long been known 
in theory |45| and also verified in experiment [46[ that 
in a one-dimensional random walk — so-called "single-file 
diffusion" — the mean-squared displacement of the parti- 
cles grows with precisely the same exponent as in our 
case. However, it remains an open question whether 
in QA systems the main effect in single-file diffusion — 
nonovertaking particles — is prevalent, or if other (possi- 
bly compensating) reasons lead to a coincidental agree- 
ment in z. 



F. Trapping 

In this subsection we will elucidate some aspects of 
the geometrical structure imposed by the quenched ma- 
trix upon the fluid immersed therein. The key concept 
involved is that of "voids," i.e., domains of space that 
may be fully explored by a fluid particle if placed within. 
In HS-QA systems distinct voids are separated by infi- 
nite potential walls, and there exist voids of finite vol- 
ume. Particles located in such finite void cannot travel 
infinitely far away from their initial location; such par- 
ticles will henceforth be denoted as "trapped particles," 
and the corresponding void as a "trap." Due to the statis- 



tical nature of the matrix structure, at any nonzero ma- 
trix density there exist traps. Likewise, for any nonzero 
density of traps there will be trapped particles since the 
initial locations of the fluid particles are randomly dis- 
tributed. 

Two types of nontrivial questions remain to be as- 
sessed. First of all, there exists the aforementioned dis- 
tinct (f>^ above which all fluid particles are trapped. Dif- 
ferently phrased (cf. Sec. lIIIE| . upon varying (/),„, a per- 
colation transition of the space accessible to the fluid par- 
ticles takes place: for (pn-i < 0m) there exists an infinitely- 
large void, whereas for 0„i > 0,*^ ^^Y void is of finite size. 
Naturally, the question arises what the precise value of 
(j)^ is. In the following we will investigate the behavior 
of fluid particles that explore the voids by studying the 
distribution of their displacements, as was done before 
for other types of confinement [2l| . 

Our quantity of choice for the desired analysis is the 
self-part of the direction-averaged van Hove correlation 
function 



Gs{r,t) = {6 (r - Arit))) 



(7) 



where Ar{t) — |r*(i) — r*(0)| and by construction 
L Gs{r,t)dr = 1. Strictly speaking, to assess the void 
structure only knowledge about the infinite-time limit 
Gs(r, i— )-(X)) is required. However, since data are readily 
available, in the following we will also discuss some fea- 
tures of the time evolution of Gg (r, t) in order to corrob- 
orate and extend our findings from the previous sections. 
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Using the self-part of the van Hove function a number 
of quantities of interest can be assessed at least qualita- 
tively. For instance, the matrix density at percolation, 
(/)j^, can be estimated by exploiting the fact that as the 
system control parameter (here 0ni) is varied the average 
size of the "aggregates" (in this case the traps) diverges at 
the transition. For sufficiently small r, the distribution of 
trap sizes is reflected in Gs{r^ i— J-oo) since in the infinite- 
time limit the nontrapped fluid particles have diffused 
far away. Concentrating on a particular distance f ~ cr, 
the function Gs(r, t^oo; (f>^) will exhibit two maxima as 
(/)m is varied, one for (/)„! < (p*-^ and one for 0,„ > 0Jjj. 
The corresponding integral L Gs{r,t-^oo) dr represents 
an estimate for the fraction of fiuid particles located in 
traps. Of course, in simulations t^oo has to be approx- 
imated by some finite time t greater than the structural 
relaxation time of the system. 

In Fig . [T^ and [T51 the van Hove function is displayed at 
various f as a function of r. For visual guidance, curves 
are highlighted for selected times. The times of successive 
light (gray) lines differ by ^ VlO, i.e., there are about six 
curves per time decade. In Fig.fT^a) Gs(r, t) is shown for 
(pm = 0.05 and 0f = 0.50, a bulklike state point like those 
discussed in Sec. IIIICI At short times Gs(r, i) exhibits 
a single, approximately Gaussian peak. As t increases, 
the maximum value of this peak decreases, but notably 
its r position remains nearly constant even for large t. 
This is another manifestation of the cage effect prevalent 
in glassy systems j41| . Further support for this picture 
arises from the second peak in Gg (r, t) that emerges at 
r ~ 1.1 for t > 10^, which indicates the presence of hop- 
ping processes. For very long times even a moderate third 
peak appears at r ~ 2 as Gg {r, t) recovers a Gaussian-like 
shape centered at r > 3. 

The validity of these arguments becomes clear when 
considering Fig. \V2\h) which depicts Gs{r,t) for (f>„^ = 
0.20 and (j)i = 0.22. For this state point, represent- 
ing intermediate matrix density, similar features can be 
observed for short and intermediate times, only with 
the first peak decreasing slower in height, the second 
peak being less pronounced, and the third peak missing. 
However, for large t and small r the shape of Gs(r, t) 
is markedly different from that in Fig. [T2ja): even for 
i = 7.1 X 10"*^ (the largest time considered) we find Gs{r,i) 
to be nonzero for small ?■. Since curves of Gs{r,t) very 
much resemble each other for t < t and r < 0.5, it is rea- 
sonable to assume that in this range Gs{r, t) represents a 
reliable approximation to Gs(?', i— >oo). Hence we can in- 
fer the presence of trapped particles, with their fraction 
being x ~ J^ Gs(r, i) dr = 8% using r = 0.5. Moreover, 
since Gs{r,i) is close to zero for r > f, there likely exist 
only few traps with a spatial extent exceeding r. 

The trend seen from Fig. [T^a) to [T^Jb) is contin- 
ued by the systems represented in Fig. [T31 Panel (c) 
of that figure shows the state point at 0m = 0.25 and 
(j)i ~ 0.10, that is, close to the diffusion-localization tran- 
sition (discussed in Sec. IIIIEP : the state point in panel 
(d) is situated well in the localized phase at 0,„ = 0.30 




Figure 14: (Color online) Total intermediate scattering func- 
tion F{k, t) (lines) for selected (j}{ values at (j>m = 0.05 and 
k = 7.0. Symbols denote the right-hand side of Eq. (|10|) . 
Error bars: see Fig. [S] 



and (j){ = 0.05. The change in shape over time becomes 
less pronounced as (j)„i increases, turning the second peak 
into a mere shoulder. 

The largest times depicted in Fig. [T51 again represent a 
sound approximation of Gs(r, t—>-oo). For the state point 
in panel (d) we numerically obtain a; ~ J Gg (r, t) dr = 
98% using f = 3 and i = 7.8x10^, which is in excellent 
agreement with the figure x = 100% that would be ex- 
pected in an infinitely-large system residing beyond the 
percolation transition. For the state point in (c), using 
the same f and t we obtain x ss 58%, which, however, 
is rather imprecise since Gs(r, i) is considerably greater 
than zero. 

In Fig. [T^ c) and (d) the correlator Gs{r,t) differs 
markedly from that in Fig. [T^a) and (b). For one, in 
both (c) and (d) a second peak is present even for i — >■ oo, 
suggesting that on average traps are large enough to allow 
for hopping. Another feature is less obvious but more im- 
portant to the percolation transition: for large distances 
Gs(r, t— >oo) is greater in (c) than in (d) for the same r. 
This implies that the mean trap size is larger in (c) than 
in (d) , which is not unexpected — in fact this confirms the 
transition to take place between (b) and (d). Thus, us- 
ing Gs{r,t) we find 0.22 < 0Jjj < 0.30 as an estimate for 
the matrix packing fraction at the percolation transition, 
which is in accordance with our findings in the previous 
sections. 



IV. DISCUSSION 

In our presentation of the global aspects of the sys- 
tem's dynamics, as well as of its behavior along selected 
paths in the {(/)m,</>f} plane, we have focused our atten- 
tion on self and connected density correlators as indica- 
tors of the single- particle and collective dynamics. This 
choice is motivated by the key role played by these dy- 
namic correlation functions within the MCT framework 
for confined fiuids, which provides detailed predictions 
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Figure 15: (Color online) self-intermediate scattering func- 
tions for a series of 0m values at fixed (jii = 0.10 and k = 7.0. 
(a) self-correlator Fs{k,t), (b) self-correlator without long- 
time plateau Fs{k,t). For the two largest values of 0m in 
panel (b) it was not possible to reliably identify the long-time 
plateau. Therefore the corresponding self-correlators are left 
unshifted. 
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Figure 16: (Color online) Kinetic diagram based on F{k,t). 
Symbols: time tt needed for 7^(^=7.0, t) to decay below F* — 
0.1. Thick solid (red) line: interpolation through points for 
which it ^ 10°. 



for Fs(k, t) and Fc{k, t). Before discussing in more detail 
the comparison between the theoretical scenario of MCT 
and our simulations, it is instructive to study correla- 
tion functions related to Fs{k,t) and Ff.{k,t). Namely, 
we will investigate here the total intermediate scattering 
function 



^ 




Figure 17: (Color online) Kinetic diagram based on Fs{k,t). 
Symbols: time fs needed for Fs(fc=7.0, i) to decay below 
F* = 0.1. Dashed (blue) line: interpolation through points 
for which ie — 10^. Solid (blue) line: interpolation through 
points for which ta ^ 10^ (from Fig. [3|. Dash-dotted (red) 



line: interpolation through points for which ic 
Fig.ig). 



10" (fro 



F{k,t) 



(pk(i)P-k(0)> 



Sik) 



(8) 



as well as a modification of the self-intermediate scatter- 
ing function in which the long-time plateau (as observ- 
able, for instance, along paths II and III) is subtracted 
out. This will allow to clarify aspects related to the ki- 
netic diagram of the system and to facilitate the interpre- 
tation of the numerical results in the light of the MCT 
predictions. 

By construction, fluids confined in porous media pos- 
sess nonzero average density fluctuations, (pu) ^ 0, in- 
duced by the matrix structure. This purely-static back- 
ground becomes visible in the long-time limit of the to- 
tal intermediate scattering function F{k,t). Collective 
relaxation phenomena in confined fluids are therefore 
more conveniently described by the connected correlator, 
Eq. ([S]), in which only the fluctuations of the microscopic 
density are considered. Nevertheless it is interesting to 
also inspect the shape of F{k, t) as a function of the state 
parameters. In Fig. [T3]we show F{k,t) for some selected 
densities along path I (</)„! = 0.05) for the same value 
of k considered in Fig. El Except for the lowest fluid 
density, (j)f = 0.05, the F{k,t) correlator attains a fi- 
nite, sizeable plateau at long times. To check whether 
this long-time plateau is due to the purely-static back- 
ground mentioned above or actually reflects a nonergodic 
behavior of the system, we also include in this figure the 
function [Sc{k)Fc{k,t) + Sb{k)]/S{k), where 



Shik) = (pk> (p-k) 



(9) 
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is the so-called blocked structure factor, and Sc{k) = 
S{k) — S'b(fc) is the connected structure factor. It is easy 
to show that, if the system is ergodic, the blocked struc- 
ture factor is the infinite-time limit of the un-normalized 
total intermediate scattering function [26[ and that 



F{k,t) 



5c(fc)fc(fc,t)+5b(fc) 

Sik) 



(10) 



holds as an equality. 

From Fig. [13] we see that this equality is indeed fulfilled 
for 0f < 0.48, while slight deviations occur at the largest 
fiuid packing fractions. Note, however, that the calcula- 
tion of S'b(fc) is also affected by statistical noise (see error 
bars in Fig. [H)) and systematic effects [4J]. Nonetheless 
Eq. (jlOp provides an effective criterion to discriminate 
nonequilibrated samples as far as the collective relaxation 
is concerned. We also note that if we employ the crite- 
rion F{k=7.0,tt) = 0.1 to define the relaxation time it, 
the latter strongly overestimates the slowness of the col- 
lective relaxations compared to tc due to the emergence 
of the finite plateau in F(k,t) upon increasing (f>{. This 
should be born in mind when discussing decoupling phe- 
nomena between single-particle and collective properties. 

A similar, yet physically different long-time plateau 
is observed in the self-intermediate scattering functions 
along path II and III (see Fig. [5] and [TU)) . As discussed 
in Sec. IIIIFl the finite value of Fs(fc,i— >oo) reflects the 
trapping of particles in disconnected voids of the ma- 
trix structure. Since the existence of a rising plateau 
affects the evaluation of the relaxation times tg, it is sen- 
sible to define a modified correlator Fs{k,t) = Fs{k,t) — 
Fs(fc,i— >-oo), where the long-time limit Fs{k,t^(X)) is 
identified by some construction [51| . The results of this 
procedure are collected in Fig. [12] for selected state points 
at (pf = 0.10 (path III) and compared to the correspond- 
ing Fs{k,t). At the two largest (pm the relaxation to the 
finite plateau is not complete within the time range con- 
sidered, and the correlators are left unshifted. Evalua- 
tion of the single-particle relaxation time ig, defined by 
Fs{k=7.0,ts) = 0.1 as for the other correlators, allows to 
filter out the contribution to the single-particle dynamics 
due to trapping. 

We constructed two alternative kinetic diagrams for 
the system at hand using the relaxation times it and 
is obtained from F{k,t) and Fs{k,t). As in Sec. IIIIBl 
lines drawn in Fig. [TBI and [TTl correspond to a fixed, con- 
stant value of the relaxation times (isorelaxation times 
lines). Let us first focus on Fig. [161 displaying the ki- 
netic diagram based on it. In this case, the range of state 
points available to determine the isorelaxation times lines 
is more limited due to the incipient growth of the plateau 
in the region of large (pi and large (pm ■ Nevertheless it is 
clear that the shape of the estimated arrest line is non-re- 
entrant, confirming the analysis based on the connected 
correlators (see Fig. [S]). 

A more interesting effect is observed in Fig. [T7| where 
we show the kinetic diagram based on ig. For better 



comparison with the analysis in Sec. IIIIBl we also in- 
clude the isorelaxation times lines obtained from the self 
and the connected correlators. Interestingly, the shape 
of the kinetic diagram for ig strongly resembles that ob- 
tained from the connected correlator. In particular, the 
iso-ig lines do not bend rapidly toward (p{ = as (p^ 
is increased but rather stretch to larger (p^ values, fol- 
lowing the trend of the iso-ic line. This effect can be 
understood as follows: after subtracting the long-time 
plateau, which is due to trapping, the residual relaxation 
of Fs{k, i) is governed by a weak, collective caging effect. 
The relaxation of Fs{k,t) toward the long-time limit oc- 
curs, in fact, on a similar time scale as the relaxation of 
Fc(fc, i). This supports our interpretation of the complex 
relaxation pattern of Fs{k,t) for large matrix packing 
fractions (see Fig. [5]) as a superposition of trapping and 
caging effects, and confirms the connection between the 
diffusion-localization transition and the appearance of a 
nontrivial decoupling between single-particle and collec- 
tive properties |3l|. 



V. CONCLUSION 

In this contribution we have investigated the dynamic 
properties of a fluid conflned in a disordered porous ma- 
trix. In our investigations we have reduced such system 
to a simple QA model where the matrix is a quenched 
configuration of a liquid and the fiuid particles equili- 
brate within that matrix. Both matrix and fluid parti- 
cles are chosen to be hard spheres of equal size. Thus, the 
parameter space characterizing the system reduces to a 
two-dimensional plane spanned by the packing fraction of 
the matrix particles and that of the fluid particles. Inves- 
tigations have been carried out via event-driven molec- 
ular dynamics simulations, observables were evaluated 
using the double-averaging recipe characteristic for QA 
systems. Specifically, we have evaluated the connected 
and the single-particle intermediate scattering functions, 
the mean-squared displacement, and the self-part of the 
van Hove function. 

The data and conclusions summarized here repre- 
sent a counterpart to the recently presented results ob- 
tained from a theoretical framework J25l - l27t . The latter 
approach — an extension of MCT to the QA protocol — 
had predicted a number of intriguing features in the ki- 
netic diagram. Our investigations have complemented 
this scenario of ideal transitions, unveiling the proper- 
ties of the real transitions: (i) At low matrix pack- 
ing fractions we have undoubtedly confirmed the pre- 
dicted type-B glass transition scenario, (ii) For inter- 
mediate matrix packing fractions we have found at least 
two arrest mechanisms of different nature, with the self- 
intermediate scattering function displaying a more com- 
plex behavior than its connected counterpart. This sce- 
nario can be rationalized as a superpositioning of the 
effect of trapping in disconnected voids and of a weak 
caging mechanism collective in nature, (iii) Varying at 
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Figure 18: (Color online) Solid (red) line: probability 
P{N, 6) = 50% to set up a system instance with A^ = 
Nm+N{ ~ 1000 particles within 6* = 10 CPU minutes us- 
ing the custom algorithm described in Appendix [X] Lower 
left: PiN,9) > 50%, upper right: P{N,e) < 50%. Dotted 



(green) line: 
and [111). 



Sr^{t) dynamic arrest criterion (see Sec. IIIIBI 



low fluid packing fractions the density of the matrix we 
found evidence for a diffusion-localization transition, re- 
flected in the self-intermediate scattering function as well 
as in the mean-squared displacement and in the self-part 
of the van Hove function, (iv) Despite considerable ef- 
fort we were unable to support the predicted re-entrant 
transition scenario. 

In order to obtain deeper insights into the open ques- 
tions we have additionally considered along path II the 
total intermediate scattering function, and along path III 
we have subtracted from the self-intermediate scattering 
functions the long-time plateau as there is evidence that 
the latter is due to trapping. Our analysis of these data 
supports the previously established notion of an interplay 
between the two phenomena dominating the arrest tran- 
sition at intermediate matrix packing fractions, namely 
trapping and caging. Work to disentangle the effects of 
the latter on the dynamic correlation functions is cur- 
rently underway. 
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Appendix A: System setup details 

As mentioned in Sec. |TT] the initial setup of a HS-QA 
system represents a formidable challenge. In this ap- 
pendix we will point out algorithms known to succeed in 
setting up disordered one-component HS systems at high 
density, comment on problems in extending those algo- 
rithms to the QA protocol, and report on our custom 
solution to the problem. 

The challenge is the following: Given the positions of 
the A^ni immobile obstacle particles, permissible positions 
are sought for the A^f particles that during MD will be 
allowed to move. Unfortunately, to this respect multiple 
complications arise in HS-QA systems, the three most 
important of which are: (1) particle overlaps are strictly 
forbidden, (2) in the resulting configuration the matrix 
particles have to occupy precisely-specified locations, and 
(3) all voids — notably also the disconnected ones — have 
to be considered as locations for the fiuid particles. 

The simplest method — trial-and-error insertion of the 
fiuid particles — is useless already for setting up bulk high- 
density systems of hard spheres. For the latter task a use- 
ful straight-forward method (although prone to crystal- 
lization) is compressing a low-density system, for instance 
by coordinate re-scaling combined with random particle 
displacements, or by applying a unidirectional force in 
a bounded system. However, there also exist elaborate 
algorithms to maximize the density in disordered HS sys- 
tems; they employ for instance serial deposition [37[ or 
overlap elimination |38| . 

Unfortunately, the complications mentioned before 
render these methods difficult to adapt and/or ineffec- 
tive. Therefore, in order to efficiently achieve high den- 
sities in HS-QA systems we devised a custom algorithm, 
which consists of the following three steps: 

(1) Upon insertion of a fiuid particle its "real" diame- 
ter (the diameter to be used during MD simulations) is 
reduced to some minute value ("defiated") so that the 
probability of an insertion via trial-and-error is greatly 
increased. This way, throughout the search of an allowed 
configuration all particles can be present in the system 
simultaneously. 

(2) In order to find such allowed configuration a sim- 
ple Metropolis Monte-Carlo algorithm [43] is employed 
using the HS potential. Since using this algorithm with 
the minute diameters is meaningless, at the beginning of 
each sweep (one displacement attempt for each fluid par- 
ticle) the diameter of each deflated particle is increased 
( "inflated" ) to the maximum value possible without over- 
lapping another particle. If for a particle the latter value 
is greater than its real diameter then inflating is done for 
this particular particle and it is assigned its real diameter 
("fully inflated"). 

(3) A disconnected void may be filled with (almost) 
any number of deflated particles; however, the particular 
void may be too small to accommodate all of these parti- 
cles if they were fully inflated. To remedy this problem, 
before certain sweeps all particles that are not yet fully 
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inflated are removed from the simulation box and then re- 
deflated and re-inserted according to step 1 . The number 
of sweeps separating two such removal-and-reinsertion 
procedures is gradually increased so as to allow for any 
number of sweeps to fully inflate a particle while still 
quickly "draining" crammed voids. 

Obviously, as soon as all particles are fully inflated 
the setup routine is completed and a configuration fulfill- 
ing the QA requirements has been found. In step 2 the 
fluid particles arc inflated serially. This leads to the dis- 
tribution of diameters covering a wide range at all times 
during iteration, except when close to the approach of an 
allowed configuration. Since this introduces a "fluctuat- 
ing polydispersity," we found the generation of ordered 
states at high densities to be strongly suppressed even 
for bulk monodispcrse systems. 

Figure [18] gives an overview over the state points at 
which, according to the above-described steps, a system 
instance with N = Nm+N[ ~ 1000 particles may be set 
up with probability P(iV, 0) ~ 50% within 9 = 10 CPU 
minutes (Intel Core 2 Duo E8400). Comparing the so- 
found line of discrimination (red solid line in Fig. llSp with 
the Sr^{t) dynamic arrest criterion (green dotted line in 
Fig. mi see Sec. IIII Bl and [Slf ) and other features of the 
kinetic diagrams, one can see that even using our custom 
routine the system setup represents a challenge to access 
the state points of interest. 

Note that in practice a probability of 50% is unsatis- 
factory. Ideally, none of the setup attempts should be 
rejected since as little bias as possible should be exerted 
on the statistics. However, already at relatively low (/)„! 
achieving P = 100% is strictly impossible, the reason be- 
ing that there exist peculiar matrix configurations that 
prohibit all fiuid particles to be inserted. Nonetheless, by 
prolonging the CPU time 9 it is possible to increase P 
since more theoretically possible system configurations 
can actually be realized. Throughout this paper only 
state points have been considered for which P > 90% 
was achieved. 



Appendix B: Finite-size effects 

As mentioned in Sec.|lTl throughout this work we chose 
to assign the same constant total number of particles N = 
N„^+Ni = 1000 to most systems simulated. The purpose 
of this appendix is to determine for which state points 
this number is not sufficient to rule out the presence of 
finite-size effects, and to describe our pertinent solution. 

Strictly adhering to the approach of constant N, state 
points with a large disparity of 0m and (j)[ attain low val- 
ues of either TVm or N{. Neither is desirable since it will 
render the calculation of properties pertaining to either 
the matrix or the fiuid inaccurate. In this work, the most 
extreme case for low N„^ is the state point {cpni = 0.050, 
(j)f = 0.505), with N = 1000 entailing N^ = 90. The 
corresponding extremal case for N{ is the state point 
{4'ni = 0.3000, (j)f = 0.0025), leaving only eight fluid par- 




Figure 19: (Color online) Assessment of the finite-size effects 
present at (pm = 0.2500 and (j>f = 0.0125. (a) Fc{k,t) and (b) 
Fs(/c, t), both displayed for k € {1, 2, 4, 8} as well as iV = 1000 
(red solid line, squares), N = 4000 (green dashed line, circles) 
and N = 10 000 (blue dotted line, triangles). 



tides to move in a matrix of 992 particles. Hence, in the 
present work the much more severe problem lies in low 
values of A^f . 

To remedy this problem we chose to adjust A^ such that 
neither N^ nor A^f assume a value below a lower limit of 
A^* = 50. The previous brief analysis shows that in none 
of our systems A^m < A^*, which Unfits the discussion 
to state points with low A^f. In the following, we assess 
whether or not the choice of A^* is reasonable. 

In order to limit the effort, we consider a single state 
point representative of the problem of low A^f. Specifi- 
cally, we inspect various observables at 0m — 0.2500 and 
(j)f = 0.0125, i.e., at a state point close to both (j){ ~ 
and the percolation transition in the voids (cf. Sec. IIII El 
and lllTF)) . Using A^ = 1000 this system attains A^f = 48 
which is reasonably close to A^* . To investigate finite-size 
effects we computed Fc{k,t), Fs{k,t), Sr'^{t), and z{t) for 
A^ = 1000; 4000; and 10 000 while in each case averaging 
over ten matrix configurations. 

In Fig. [Tithe dependence of Fc{k, t; A^) and F^^k, t; N) 
upon A^ is visualized. The intermediate scattering func- 
tions arc shown for selected wave vectors, where any 
finite-size effects should be most prominent at small k. 
Indeed; for fc — 1.0 and t > 10^ there are noticeable dis- 
crepancies between Fc(fc, t; 7V=1000), Fc(fc,i; 7V=4000), 
and Fc(fc, t; A^=10 000); however, this is hardly surpris- 
ing considering that even for A^ = 10 000 the corre- 
sponding length scale 27r/fc = 2tt is only slightly smaller 
than half the edge length of the simulation cell. In 
fact; already for k — 2, there is no distinguishable dif- 
ference between Fc{k,t;N ^1000), F^k, t;N =4000), and 
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Figure 20: (Color online) Assessment of the finite-size effects 
present at 0^ = 0.2500 and (j)i = 0.0125. (a) cSr^(i) and (b) 
z{t), both displayed for A'' — 1000 (red solid line, squares), 
iV = 4000 (green dashed line, circles) and iV = 10 000 (blue 
dotted line, triangles). 



Fc(A;,i;^=10 000). On the other hand, Fs{k,t;N) dis- 
plays small consistent discrepancies between N = 1000, 



4000, and 10 000; however, the worst deviation between 
Fs{k,t; iV=1000) and Fs (fc, i; iV=10 000) (which is at fc = 
8.0 and t = iniax = 2.25x10^) is no more than ^12%. 
It is not surprising that the Fs{k,t; N) curves in this 
case exhibit larger deviations from one another than the 
Fc{k,t;N) curves considering that the latter approach 
zero rather quickly whereas Fs{k,t; N) does not com- 
pletely relax. Also, we point out that any discrepancies 
in Fs{k,t;N) for different N values are well within the 
range of the error bars. 



Similarly, Fig. 



illustrates the dependence of 



Sr'^{t;N) and z{t;N) upon N. The mean-squared dis- 
placement displays only minor differences for N = 1000, 
4000, and 10 000: The discrepancy between Sr^{t; 1000) 
and fr^(t; 10 000) at t = imax amounts to ~24%; how- 
ever, since throughout this work we are interested in 
Sr'^{t;N) on logarithmic scales, the above deviation can 
be considered sufficiently accurate. z{t; N), on the other 
hand, is clearly more susceptible to stochastic errors than 
(5r^(i; N), leading those errors to dominate finite-size ef- 
fects in z{t>3.6xl0^; N=1000). For t = 3.6x10^, we 
find a difference of mere 17% between z{t; A^=1000) and 
z{t; A^^IOOOO), which permits almost quantitative inter- 
pretation. 

We conclude that A^f = N* is sufficient to yield reliable 
results for all observables considered. None of the latter 
exhibited a deviation of more than 20% between N = 
1000 and the tenfold value for N, which validates all 
conclusions presented in this work. 
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More detailed, the procedure to set up a QA matrix reads 
as follows: (1) define an interaction between the matrix 
particles, (2) find some allowed configuration for a system 
containing just the matrix particles, (3) equilibrate that 
system, (4) take a snapshot at some instant of time, (5) 
use the particle positions therein as obstacle positions. 
Increases or decreases of A'^^ and/or A^f to accurately 
approximate the desired fraction N^n/Nf resulted in slight 
deviations from the standard N = 1000. 
The denominator, although equal to unity, has been re- 
tained to unveil the structure of the correlator. 
The numerical procedure to evaluate Fs{k, t) for fixed k is 
as follows: first, the global minimum M = mint>o Fs{k, t) 
of the self-intermediate scattering function is identified. 
Then, the standard deviation, S, of Fs{k,t) is computed 
over different time ranges {imin, imax}, keeping the upper 
bound imax of the range fixed at the maximum time avail- 
able while progressively increasing the lower bound tmin. 
If S/{1 — AI) is below some preset tolerance threshold then 
Fs{k,t) is considered to be constant in the corresponding 
time range. If this condition is met while tmax/imin > 10 
then Fs{k,t) = (Fa(fc,i)-M)/(1-M). Otherwise, Fs{k,t) 
has not reached its long-time plateau within the time 
range considered. 



